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LEGENDRE  METHODS  ON  CHEBYSHEV  POINTS 


Wai  Sun  Don  and  David  Gottlieb^ 
Division  of  Applied  Mathematics 
Brown  University 
Providence,  RI  02912 


ABSTRACT 

We  present  a  new  collocation  method  for  the  numerical  solution  of  partial  differential 
equations.  This  method  uses  the  Chebyshev  collocation  points,  but  because  of  the  way  the 
boundary  conditions  are  implemented,  has  all  the  advantages  of  the  Legendre  methods.  In 
particular  Li  estimates  can  be  obtained  easily  for  hyperbolic  and  parabolic  problems. 
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1  Introduction 


Polynomial  pseudospectral  (or  collocation)  methods  have  been  extensively  used  in  the  nu¬ 
merical  solutions  of  partial  differential  equations.  The  underlying  idea  in  those  methods  is 
to  approximate  the  unknown  function  by  an  interpolation  polynomial  at  some  pre-described 
(collocation)  points.  The  polynomial  is  then  required  to  satisfy  the  PDE  at  the  collocation 
points.  This  procedure  yields  a  system  of  ordinary  differential  equation  to  be  solved. 

Historically,  (see  [10])  the  first  such  points  to  be  used  were  the  Chebyshev  collocation 
points 

Xj  =  cos(^)  Q<j<N. 

Those  points  were  chosen  because  they  allowed  the  use  of  Fast- Fourier- Transforms  in  the 
computatiouo.  U  was  only  later  (see  [7])  that  those  points  were  idefitified  with  the  nodes 
c'f  thv  Gauss- Lobat to- (G-L-C)  quadrature  formula.  This  observation  is  the  key 
In  the  stability  analysis  of  the  pseudospectral  Chebyshev  methods.  The  G-L-C  quadrature 
formula  led  to  the  weighted  La  norm 


However,  it  has  been  noted  in  [8]  that  this  is  not  a  natural  norm  for  hyperbolic  equations. 
In  fact  the  differential  equation  is  not  well  posed  in  this  norm.  Also  it  complicat('*d  the 
stability  analysis  even  for  parabolic  equations.  The  theory  (and  therefore  the  confidence  in 
applying  those  methods)  is  not  complete. 

Once  the  connection  between  the  collocation  points  and  the  Gauss  Lobatto  points  is 
established,  it  is  natural  to  use  die  nodes  of  the  Gauss- Lobatto- Lcjfrndrr  (G-L-L)  quadrature 
formula.  We  refer  the  reader  to  [2]  for  review  of  those  methods.  Recently  [Ij  an  0(N  log  N) 
method  was  proposed  for  the  Legendre  points.  The  main  problein  with  those  points  are  that 
they  are  not  given  explicitly,  and  their  evaluation  for  large  AT  is  not  iQobust  due  to  roundoff 
errors. 

In  this  paper  we  present  a  method  (and  name  it  Thr  ('lirbytJttV'Lrgrndrr  Mrthod)  that 
ha.s  the  advantages  of  both  the  (Chebyshev  and  Legendre  methods.  The  method  utilizes  the 
('hehyshev  collocation  points  allowing  the  use  of  fast  Fourier  algorithms  and  avoiding  the 
roundoff  error  associated  with  computing  the  Legertdre  grid  points.  The  boundary  conditions 
nre  imposed  via  a  new  (mnalty  technique  in  such  a  way  that  the  method  is  stable  in  the  usual 
Li  norm  (ratlrer  than  the  weighted  Lj  norm).  Hen*‘e  the  (.‘hehyshev- Li^^endre  method  enjoys 
the  advantages  of  the  Chebyshev  method  as  well  as  those  of  the  Legendre  tnethod. 


The  implementation  of  the  boundary  conditions  is  done  by  a  penalty  method.  A  penalty 
term  is  added  to  the  F’DE  at  all  grid  points  in  such  a  way  that,  in  the  limit  of  number  of 
grid  points  tend  to  infinity,  the  boundary  conditions  are  satisfied.  This  procedtire  seems  to 
be  better  than  the  direct  imposition  of  the  boundary  conditions,  and  in  our  case  has  the 
extra  advantage  of  yielding  the  Legendre  method  at  the  ('hebyshev  points. 

A  similar  idea  had  been  tried  by  Reyna  (11).  The  difference  between  his  approach  and 
ours  is  in  the  imposition  of  the  boundary  conditions.  Instead  of  transforming  from  the 
(/’hebyshev  basis  to  the  Legendre  one  as  in  [11],  we  impose  the  boundary  conditions  via 
penalty  method  and  through  that  automatically  switch  to  the  Legendre  basis  mthout  tising 
it  in  the  diffenntiation  pit>eedurr. 

The  paper  is  organized  as  follows: 

In  .Section  2  we  quote  the  essential  formulas  for  the  use  of  (Ihebyshev  and  Legendre 
methods. 

In  Section  3  we  present  the  (’hebyshev-Legendre  method  for  hyperbolic  equations.  In 
subsection  3.1  we  describe  the  method  and  prove  (in  Theorem  (3.1.1))  an  energy  estimate 
to  show  stability.  In  Theorem  (3.1.2)  we  bring  another  version  of  this  energy  estimate. 

In  subsection  3.2  we  consider  the  relationship  of  the  new  method  to  the  Legendre  penalty 
method  and  show  that  the  differentiation  matrices  of  the  two  methods  are  related  via  a 
similarity  transformation.  This  fact  is  proved  in  Theorem  3.2.2 

In  Section  4  we  discusis  and  prove  the  stability  of  Chebyshev- Legendre  methoii  for  the 
heat  equation  with  Robin  boundary  conditions. 

Section  5  concludes  the  paper  with  some  numerical  experimentations  with  the  new 
method. 

In  future  work  we  will  report  on  the  convergence  results  of  the  new  method  for  nonlinear 
hyperbolic  equations. 

2  Preliminaries 

This  Section  is  devoted  to  the  definitions  of  the  pseudospectra)  methods  to  be  used  later.  We 
will  discuss  (Hiebyshev  and  Legendre  methods  which  are  based  on  the  (^rbyskev  polynomiak 

r/tf(jr)  =  cos(A'cos"*  x)  (2.1) 

and  the  Lrgmdrr  polynomuds 

respectively.  .Associated  with  these  two  polynomials  are  .several  tiauss-lype  quadrature  for> 
mulas.  We  will  consider,  in  this  paper,  the  Canss  iohotlo  type  formulas. 
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We  start  by  defining  the  Chebyshev  Collocation  points  by 

Xj  =  cos(^)  0  <  j  <  (2.3) 

These  points  are  the  zeroes  of  the  polynomial  (1  —  x^)Tj^{x)  and  associated  with  it,  we  have 

The  Gauss  Lobatto  Chebyshev  Quadrature  Formula; 

Let  f{x)  be  a  polynomial  of  degree  2N  -  I,  then 


}=0 


where  the  weight  Cj  are  given  by 


c,  =  ^  l<j<iV-l 


Co  —  Cs  — 


Similarly,  the  Legtndrr  Collocation  points  yj  are  defined  as  the  roots  of  the  polynomial 
(1  —  x'*)/^(x).  For  these  points  we  have 

The  Gauss  Lobatto  Legendre  Quadrature  Formula: 


Let  /(x)  be  a  polynomial  of  degree  2iV  -  1,  then 


EasijK"  /'/(«)<« 


where  the  Gauss  Lobatto  weights  >jJj  arc  given  by 


mA)  —  W'.v  — 


;V(.V  +  1) 


t'nlike  the  Chebyshev  points  that  are  known  explicitly,  there  is  no  explicit  forntula  for  the 
Legendre  points  y^,  they  have  to  l>e  computed  numerically.  It  is  int«*resting  though  that  iltere 
is  a  sintpie  formula,  easily  and  robustly  computed,  for  the  values  of  the  Legendre  polynomials 
at»d  their  derivative  at  the  Chcbyshtc  points.  In  fact  we  have  the  following  explicit  forntula 
for  /xUj),  (taken  from  (4).  page  180) 


>‘-t 


P'^4cos0)  ss  U*?  ^  7  ~  “•  I  “ 

^  m!(.V  -  I  -  tn)! 


•  I 


Pseudaspectral  (or  (V)llocatiou)  methods  are  based  on  interpolations  at  the  points  Xj  t)r 
Pj.  (/'onsider  the  polynomials 

QU^)  = 

Qc(t)  =  (l-x’)n(x), 

and  define  the  Legendre- Lagrange  polynomials  by 


h,ix)  = 


Ql{x) 

(•r-tfj)QL(yj)’ 


and  the  (’hebyshev- Lagrange  polynomials  by 


Qc(^) 


Then,  the  Legendre  interpolation  operator  Ii  is  definerl  by 

(/L/)(-r)  =  i;/(y,)/»j(^) 


(2.10) 


(2.11) 


whereas  the  (’hebyshev  interpolation  operator  Ic  is  definetl  by 


(/c/)(.')-EM)y,(x). 


By  definition,  we  hav« 


(2.12) 


UiDiyj)  =  /(y,) 

=  /(-r,). 

From  tlie  definition  of  the  interpolation  operators  /t  anti  Ic  we  get  the  spfrtml  dijjfrrtn- 
nation  matricta  ?>;.  and  Pc  as  follows: 


The  Pseuduspe^  tral  l.eg«*ndre  Diferentiatiun  Matri.\  P/.,  is  definetl  by 

(Pt)M  =  ^(yj). 

The  Psetidospeeiral  ( 'hebyshev  Differentiation  Matrix  Pc.  is  defined  by 

(see  |1|  for  explieil  expression.'*  for  the  malrinr.  1 


(2.l:i) 


f-i.H) 
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3  Hyperbolic  Equations 

3.1  Scalar  Hyperbolic  Equation 

In  this  Section  we  consider  the  scalar  initial-boundary  value  hyperbolic  equation 

Vf-=Ux  -  1  <  X  <  1  <  >  0  (3.1) 


with  the  initial  condition 


U{x,0)  =  fix) 


(3.2) 


and  the  boundary  condition 


V{Ut)=g{t)  . 


(3.3) 


The  Chebyshev-CoUocation  (Pseudospectral)  method  involves  seeking  an  Nth  degree 
x-polynuinial  u.v(j,<)  that  satisfies 


c^u^(j,/)  du^ixj) 


dt 


dx 


at  X  —  Xj  I  <  i  <  N 


(3.4) 


with  the  boundarv  condition 


«.v(l,0  =  «.v(.ro.t)  =  </(/) 


(3.r)) 


where  x.  are  deterinituti  in  (2.3)  (xo  =  1). 

Note  that  the  e<)uation  is  satisfietl  at  all  the  grid  points  except  at  the  l»oundary  point 
r  =  I  where  the  bonmlary  condition  is  sati.sfie<l. 

In  general,  the  term  ^m,v(x./)  b  evaluateil  at  all  the  grid  points,  with  the  us«*  of  either 
FFTor  tnairix-v«H  tor  multiplication  using  the  matrix  Pf.  tUpiation  (3.4)  is  then  a«lvanc»Hl 
at  all  the  grid  |M.iints.  The  value  of  the  M>lution  at  the  lunuxlary  is  then  updat***!  using  (3.’>). 

i»'  (.*ij  a  penalty  type  method  was  intro«lm-»Hl.  In  that  approach  we  still  usetHjuatiou 
(3.4)  for  the  mner  points  x,.  I  <  j  <  .V.  howrner  inslea*!  of  using  (3.5)  for  the  Iroundary, 
the  following  etpiatiou  i-s  satisHe<l 


du^il.t] 


iiUxix.t) 


iix 


Ui  “  f(u.v(I.O  -f^(M) 


l3.6) 


where  s  i.s  detertninetl  from  stability  considerations.  In  particular  it  l;a»l  l>et‘n  found  that 
.stability  follows  if 


Equations  (3.4)  and  (3.6)  can  be  combined  into  a  single  equation  by  noting  tliat  the 
collocation  points  xj  defined  in  (2.3)  are  the  zeroes  of  the  polynomial  (I  —  x^)7’jv(x).  Thus 
the  penalty  method  (5)  can  be  written  now  as 


dusixj,t)  dus(xj) 


di 


dx 


rax,  T 


(1  -h  Xj)Tf^(Xj) 

2Tj,{l) 


(u^(l,<)-(/(t)) 


(3.7) 


for  j  =  0, . . . ,  N. 

The  main  difference  between  the  penalty  method  (3.7)  and  the  usual  (’hebyshev  method 
given  in  (3.4)  and  (3..5)  is  that  the  numerical  solution  Un{xJ)  does  not  satisfy  the  boundary 
condition  exactly,  but  only  in  the  limit  as  Af  — ♦  oo.  The  boundary  condition  is  now  part  of 
the  equation. 

Another  penalty  method,  based  on  the  Legendre  points  y,  is  presented  in  (6).  .Similar  t«‘ 
(2.6)  we  write  this  method  as 


dusix,t) 


dt 


dx 


=«.  -  T 


(1 

2^(1) 


(u,v(Uf)-y(0) 


(3.8) 


for  j  =  0, . . . ,  N. 

The  parameter  t  is  determined  by  the  stability  requirement.  Thus  the  differential  equa¬ 
tion  is  satisfied  at  the  points  y,,  j  =  I, . . . ,  A.  At  the  boundary  Xo  =  1  «ne  uses  a  combination 
of  the  boundary  condition  and  the  differential  equation. 

An  obvious  disadvantage  of  the  method  in  (3.8)  is  that  it  utilizes  the  Legendre  points. 
However  comparing  (3.7)  and  (3.8)  shows  us  how  to  utilize  the  Legendre  penalty  metho*! 
(3.8)  at  the  (’hebyshev  points. 


The  Chebyshev«Legendre  (C»L)  Method 


Let  t*fi(x)  be  the  Legendre  polynomial  of  tlegree  jV.  In  the  (‘-L  jnetho<l  we  seek  a 
polynonnal  of  degree  N  in  x  that  satisfies 


ditKi-rfJ)  _  BusjxJ) 
m  ~  dx  ' 


. . 


(3.9) 


for  j  =  0 . .V. 

Note  that  the  pf*nal!;.  tenw  different  from  zero  for  all  the  (‘hebyshev  gri»l 

V  X..  Note  also  that  applying  (3.9)  entails  the  use  of  the  diffen'iitiation  matrix  Vr  at  the 
^  points.  In  fact,  given  U3,f(Xj.()  one  finds  the  derivative  tiase«l  on  the  (‘hebyshev 

}H»ints.  and  then  add  the  penalty  tenn  with  differrnrl  weights  at  r lyry  grid  pi»im.  The  term 
/  x(-f,)  is  evaluated  tniing  the  explicit  formula  (2.8).  This  is  *lone  onc»»  and  fur  all  ft*r  any 
grid  size 


The  surprising  fact  is  that  the  (!-L  method,  though  computed  at  the  (,’hebyshev  points, 
is  stable  in  the  usual  norm,  rather  than  the  weighted  Li  norm.  In  fart  one  can  state 

Theorem  3.1.1  ;  The  L-i  Stability  of  the  C~L  Method. 

let  u,v(  J,  /)  be  the  solution  of  (:i.9).  Let  <jJj  Im*  the  weights  of  the  (lauss  Lobatto  Legendre 
quadrature  formula  and  yj  the  nodes  of  the  same  quadrature  formula.  Let  g{t)  =  0  in  (-1.3) 
and  (3.9),  then  for 

r  >  -^  =  i\{.\  +  1) 

2uio 

the  (’-L  method  is  stable  in  the  norm.  More  specifically. 


;=0  ;=0 


f{Hl{lJ){'2uor-l)  +  ui{-\J)]dt 
Jo 


(3.10) 


Proof: 


It  follows  from  (3.9)  that 


<>u.v(x,0  dusisj)  (l+j)/^(x), 

ST-  =  -nr-  -  ■ 


(:L11) 


This  is  be<-ause  both  sides  of  (3.9)  are  polynomials  of  degree  A'  that  agree  at  .V  +  1  points, 
namely  at  the  (’hebyshev  colltMation  points  x^,0  <  j  <  .V. 

VVe  now  read  (3,11)  at  the  Lryrudrt  points  to  get 


=  2^Uv(y;.0 


>)  -  'iTuAjUxd.O 


•  • 


.Since  tin*  (»auss-Lol»atto- Legendre  quadrature  furnnda  is  exact  fur  polynomials  of  degree 
i.V  -  I  it  follows  that 

=  fiuUiMkdl 

/bO 

and  thus  the  stability  estimate  (3.10)  follows.  The  Theorem  is  thus  provett. 

Nii>te  that  unlike  the  Legendre- Penalty  method  (3.H).  in  which  one  needs  to  use  the 
Legendre  points  in  the  computations,  these  |>oiuts  do  not  appear  in  the  cuittputaliuns  in 


•  • 


the  C)-L  method.  They  are  just  introdu.^ed  for  the  sake  of  the  proof.  The  artual  computatiou.s 
are  doue  using  the  ('hebyshev  grid  points  xj. 

An  energy  estimate  based  on  the  C/hebysbev  points  xj  can  be  derived  by  using  (2.1 1}  and 
(2.12)  as  follows: 


Theorem  8.1,2  ; 

Let 

(3.12) 

^=0  <a0 

with 

N 

9j{yt)oi{ykPk . 

IhsO 

where  the  (’hebyshev- Lagrange  polynomials  gj(x)  are  defined  in  (2.10). 

Then  u«(-,<)  satisfies  the  energy  estimate 

=  ll•i«(■.0)ll'  -  /'K(1,0(W  -  I)  +  (3.1») 

Jo 


Proof ; 

Equation  (3.13)  is  really  a  restatement  of  (3.10).  .Since  Uw(x.O  b  a  polynomial  of  degree 
iV  in  X,  it  can  be  represented  exactly  by 

s 

uni-rj)  =  IcUs  =  5^  Usix(J)si(x) . 

tsQ 

Thus 

,v 

=  ^Us(-rt.t)Myj)-  (3. 14) 

(st) 

The  estimate  (.'l.|:l)  follows  from  (3.10)  upon  std»stituting  (3.14)  fur  the  values  of  UsiyjJ'h 
The  theorem  is  proven. 

The  t*.L  metlaal  can  U*  vleutnl  from  many  different  |H»int->«  of  view.  Theorem  (3.1.1) 
>h(«\v.v  that  lhi.s  metlHwl.  it*  the  constant  roeffniimt  case,  is  inpiivaletlt  l»*  ll»e  la'gemlre 
|N'nalty  tttetlual  iniro*luce«l  in  faj. 

1  hits  the  C'L  method  t$  the  realisation  of  the  Legendre  method  at  the  CTteby* 
shev  points. 

We  c|i*se  this  section  by  |H*intmS  out  that  the  estiniale  t'3. 13l  euaUles  «n«e  t**  pass  t** 
s\-sinns  easily.  ,A»  tually  il  had  l»een  »lone  in  l*»j.  Ota  lias  iti  follow  the  sftine  steps.  )t  shows 
the  melhocl  is  stable  fuf  cunstanl  cc'effirirfils  systems  of  hy{M  ;l»obc  equatiotis. 


i) 


t)  • 


•  •  •  •  •  •  •  •  •• 


3.2  The  Differentiation  Matrices 


Perhapsi  more  insij^ht  can  he  gained  if  one  compares  the  differentiation  inatrix  induced  by 
the  (%L  jnethod  Vcl  ^vitli  the  differentiation  matrices  induced  by  the  Chebyshev  penalty 
method  (3.7)  Vc,  and  the  Legendre  penalty  tnetho«i  (3.S)  t>i. 

We  start  by  noting  that  the  diff»*rentiation  matrices  Vi  ami  Pr  deHne<l  in  (2.13)  ami 
(2.14)  do  not  take  into  an  amt  any  boumlary  conditions.  The  differentiatioti  matrices  in¬ 
duced  by  (3.7).  (3.f<)  and  (3.9)  are  variations  of  the  basic  matrices  Vi  and  Pf.-,  differing  only 
in  the  method  of  imposing  boumlary  <'omiition. 

Not  -surprisingly.  Vi  and  P,-  are  similar,  after  ail  both  differentiate  (xat  tly  polytHunials 
of  degree  .V.  Since  for  these  p-ilynomials.  the  operators  li  and  Ic  are  the  same,  the  matrices 
Vi  and  Pr  represent  the  same  operation  in  a  »lifferent  basis.  This  implies  a  similarity 
relationship.  .More  spe<'itically,  th.s  relationship  can  be  written  explicitly: 

Theorem  3.2.1  ; 

Let  5  be  the  matrix  wirose  elements  Sjj,  are  given  by 

•S.i  -  Aj(-rir)  (3.1.">) 

where  hj\x)  are  the  Legem  Ire- Lagrange  polynomials  detine^l  in  (2.9).  .-\gain  Xj  are  the 
(’hebyshev  points. 

Let  T  l>e  the  matrix  wlmse  dements  7',.*  are  giveti  by 

T,.k  -  fhiyk)  (3.1b) 

where  «/;(■*“'  Phebyshev- Lagrange  p<»lynomiais  are  »lehm*tl  it»  (2.10).  j/*  an*  the  Legendre 
collot'atiou  {Joints. 

Then 


S  =  7-^  (3.17) 

and 

p,  s,9PtT  (3. IS} 

Emtd 

Since  is  a  polvtnatiial  td  df’grt’e  .V  ti  is  givicn  by 

V 


♦ 


•  i 


* 


•  ••••••• 


•  i 


Substituting  x*  and  ntaking  us«*  of  the  fact  that  </j(xt)  =  fijM,  we  get 


proving  that 


Differentiating  (3.19)  we  get 


t=a 


I  =  ST. 


<=0 

However  AJ(x)  is  itself  a  polynotnial  of  degrw  N  and  therefore  it  ean  be  expres.setl  as 


M(-f)  = 


which  leads  to 


Af  Af 

5]{-fk)  =  EH  i?,(»<)^|(ym)A«{x*) 

<s0  HtsO 

The  theorem  is  thus  proved. 

VVe  will  show  now  that  the  ditferentiatiun  matrix  induce<l  by  the  <^L  meiiuMl  Vc'i  is 
similar  to  the  ditrereutiatiuu  matrix  Pcf>  iuducetl  by  the  U'gendre  penally  metluKl  {-i.S). 
This  will  demon.strate  the  fact  that  the  (‘-t  methoil  is  the  realization  of  the  l^ginairt*  method 
on  the  (’hebyshev  grid. 


•  i 


L«*t  P(‘i  Ix'  the  differentiation  matrix  induceti  by  the  < ‘hebyshev  l,t^endrt»  methixl  (3.9) 
and  Vif  the  diR^ereutiatiot)  method  inducr.l  by  the  Legendre  penalty  ntHhiel  (:LS)  then 


Vrt  ~  SDit'T 


{im 


where  the  matrices  S.T  dtdine*!  in  (-l.l'i).  (3.16)  are  tin*  transWmatiun  matrices  l»e{wtvn 
the  ('hebyshev  points  and  the  Legendre  points. 

Note  fhai  tire  differential ioti  matrix  Pit  i* essentially  the  matrtx  D|  intrcelnip«i  ir  fi.l.il. 
modihefl  l«  take  into  ace^amt  the  (iunndary  ru»HlrtH»n»,  im|Hwril  via  jamaltj  in  {d.Xj  Thus 


(-UU 
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•  I 


In  the  same  niaiiner  we  (’.an  write  explicitly 


(T>cL)j,k  =  {'Dc)j,k  -  T 


{l+Xk)Pk{Xk) 

2/yi) 


^j.o- 


(3.22) 


Equation  (3.22)  is  a  direct  consequence  of  (3.9).  Note  that  the  full  first  column  of  Vc  is 
modified,  and  not  only  the  first  element  as  in  (3.21). 

We  proceed  by  writing  explicitly  the  elements  of  the  matrix  SVipT .  In  fact 


[SVu>T),m  -  [SVl7),m  - 


N  N 

1=0  m=0 


Thus  using  (3.18)  we  get 


[SVu>7)j.k  =  {'Pc)j,k  -  Tgj{}jo)hQ{xk) 


(3.23) 


From  (2.9) 


ho{xk)  = 


{l+Xk)PN{Xk) 

'^P‘M 


and  since  j/q  =  =  1  djiVo)  =  so  the  right  hand  side  of  (3.23)  is  exactly  the  same  as 

this  of  (3.22). 

Thus  (3.20)  is  established.  The  proof  is  completed. 


4  Parabolic  Equations 


In  this  Section  we  present  the  (Miebyshev- Legendre  method  for  the  parabolic  equation 


du 

dt  dx^ 

with  Robin  boundary  condition 


-  1  <  jr  <  1,0  0 


au{[J)  +  ihiA\,t) 

7u(-l,0  +  ^Mr(-Ut)  =</~(0- 


(4.1) 


(4.2) 


We  will  assume  that  a,  li,  7  are  non-negative  and  /» is  nun- positive.  This  a.ssun's  t  he  time 
decay  (or  non-growth)  of  u{x,t). 

We  note  that  by  now  there  is  a  very  limited  stability  tlunny  for  the  (‘hebyshev  method. 
Ill  fact  ‘  tability  had  been  proved  iirst  for  the  l)irichh‘t  ca.se  .i  =  0.  A  =  0.  (siv  [T])  and  then 
for  Neumaim  case  o  =  =  U.  [9].  Here  we  pres«'nt  the  C-L  nu'thod  ami  prove  stability  for 

the  approximation  to  (  l.l),  (4.2)  for  the  general  Robin  ca.st*. 


11 


•  ••••••• 


Denote  by  Vn  the  finite  dimensional  space  of  polynomial  of  degree  at  most  N.  We  define 
the  operator  A 

A  :  Vfi  — »  Vn 


where 


R{x,t)  =  TQQ-^{x)[B'^it)  -  (7+(J)]  +  Tf,Q~ {x)[B~ {f.)  - 


q-(x)  =  •  »■(')= 

The  numbers  Tq,  r,v  will  be  determine  later  to  assure  stability.  We  define  also  the  following 
scalar  product 

.V 

(e.  u').\  =  X!  Hyjhiyj)^'j  (•!••'') 

where  y_,, are  the  Legendre  points  anti  weights  respet'tiveiy. 

The  Chebyshev-Legendre  Method  for  Parabolic  Equations 


We  seek  the  ptdynomial  of  tlegrtv  .V  in  x.  e(x, /)  that  satisfies 


\x=,,  -  B(Xj.t)  0<j<.\ 


where  Xj  are  t'hehyshev  t'ollucation  points. 

Note  that  again,  the  work  is  tluiie  on  the  (‘fiebyshev  (loints  Xj,  the  penally  values 
are  computetl  by  ('i.?')  au«l  are  mmzero  for  a»y  r,. 

To  prove  the  stahility  of  we  s«*i  y*tf)  ~  y~{i)  -  U.  In  the  fidtowing  lemma,  we  will 
fiiul  comfit  ions  on  Tq  ami  r.v  such  that  the  opt^raiur  A  to  (»«•  semi  tH»uml»*tl. 

Lemma  4.1  ; 


12 


Let  and 


=  ^^[(l+2ic)  +  2v'/c  +  /c2] 

^((1+2k)-2V/c  +  /c«] 

with  K  ~  u>oalb. 

Let  the  operator  A  be  defined  in  (4.3).  Then 

{Av,v)n  >  ^  vi(yj)ui  (4.7) 

i=i 

provided 

T~^  <To  <T*^  (4.8) 

r-j4,  <Ts  <T+„.  (4.9) 

Proof ; 

Since  the  C>auss  Lobatto  Quadrature  formula  (2.6)  is  exact  for  polynomials  of  degree 
2N  —  1  and  since  u(i,<)  is  a  polynomial  of  degree  N,  we  have 

=  /  v(x)v„(x)dx 

=  t'(l)v«(l)  “  v(-l)Vr(~l  )  -  J  ^  v^(x)v,{x)dx 

using  the  standard  integration  by  parts  te<'bni<jue. 

Using  again  the  (Jauss  Lobatto  formula,  one  would  get 

\  N 

-  H  - »’( 1  )‘’i( I )  +  ‘d ^  1  )‘v( “  1 ) 

j=0  ;a0 

=  (^‘0) 

,.bI 

i'^(  I  )«Jta  +  t'ii  - 1  I )  ■<■  e(  - 1  )«•,(  -  1 ) 

rh»>  making  use  of  (  I.IO)  we  t  an  write 

V-j 

Ue.i-K  ==  f(l.o..Lm+ (LIU 


where 


F{x,a,b,k)  =  x{  r,,buJk  -  l)t;(a:)tu-(x)  +  TkaujkV^{x)  +  0Jkvl{x) .  (4.12) 

In  order  for  A  to  be  positive  we  need  to  choose  Tq  and  such  that  and 

F(-1,7,  1^1,  A^)  are  non-nej^ative.  For  F(l,«,/:^,0)  to  be  positive,  we  need 

(ro/itJo  -  1)^  <  4aTou;o 


oi 


TQ/i^u?o  —  2toujo{(^  +  2rtu>o)  +  1  $  0 

Tins  To  has  to  lie  between  the  roots  of  the  parabola  described  in  the  left  hand  side, 
namely  and  rl^. 

The  same  kind  of  consideration  holds  for  r.v.  Thus  F(l,a,d,0)  and  F(— 1,7,  A^)  are 

non-negative  for  the  range  of  Tq  and  given  in  (4.8)  and  (4.9),  respectively.  (4-7)  follows 
from  (4.11). 

Remarks 

1.  The  Dirichlet  boundary  condition  for  x  =  1  is  obtaineti  from  (4.2)  by  setting  n  = 
!,  li  -  0.  In  this  <•  ‘se 

r+„  =  00 

=  C?' 

which  yields  the  condition  for  the  penalty  amplitude 

lo 

2.  The  Neumann  boundary  conuition  foi  x  =  1  corresponds  to  the  case  e>  =r  0,  d  =  1.  hi 
this  case  T,t^  =  Tq”,  yielding  the  condition 

_j__N(£+_n 

2 

We  are  now  reatly  to  stale  the  stability  lh<*orer  ^or  the  ('-L  meth'id  when  applie<i  to 
paraiiolic  ^uatious  with  Hobin  bouiulary  cunditiou,s  : 

u 


t 


Let  To  and  satisfy  (4.8)  and  (4.9)  respectively.  Let  v{x,t)  €  Vs  Le  the  (J-L  approxi¬ 
mation  to  u(x,  t),  obtained  by  (4.6).  Assuming  that  </■*■(<)  =  =  0,  satisfies  the 

energy  estimate 

.tN-\ 

iv{x,t),v{x,t))s  <  {v{x,0),v{x,0))s -2  (4.13) 

j=i 

where  the  scalar  product  {f,g)N  is  defined  in  (3.5). 

Proof : 

Since  (4.6)  holds  for  j  =  G, ...,  N  and  since  v,t>xx  and  H  are  polynomials  of  degree  at  most 
yV,  we  conclude  that  both  sides  of  (4.6)  are  equal  not  only  at  the  grid  points  but  also  for 


every  x. 


dv{x,t)  d'^v{x,t 


-Rix,t)  -l<x<l 


where  R{x,  t)  is  defined  in  (4.4). 

Noting  the  definition  of  A  in  (4.3),  we  get 


Using  Lemma  4.1  yields 


=  Av. 


(v,-^)n  =  iv,Av)s 


Id  . 

2^(tM')/v  <  -  E 


and  integration  yields  the  stability  result  (4.13). 

We  stres.s  again  that  the  Legendre  coUu«'ation  points  tjj  are  "ghtist  points”,  which  are 
never  usetl  in  the  computations  but  only  in  the  priKjf  of  the  stability.  Actually  we  could 
restate  the  proof  in  terms  of  the  (3iebysbev  collocation  points  Xj  iis  in  Thetireni  3.1.2. 

5  Numerical  Results 


Case  1:  Linear  scalar  PDE 
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In  this  Section,  we  will  consider  some  numerical  examples  that  verify  our  claims  stated 
in  previous  Sections.  Consider  the  scalar  linear  initial-boundary  value  hyperbolic  PDE 


Ui  =  (4 


-I  <a;<  I,#>0 


with  initial  condition 


U{x,  0)  =  sin(2jril‘x) 


and  boundary  condition  at  x  =  I 


f/(l,<)  =  g{t)  =  sin(25rfc(l  -|- 1)) 

We  seek  an  N  degree, x-polynomial  v(x,t)  that  satisfies 

^v{Xyt)  =  Dv{Xj,t)  -  TQ{Xj){v{\,t)  -  g{t))  (5.2) 

at  (fiiebyshev  collocation  point  Xj  =  cos{irj/N),j  =  0, ...,N  and  D  is  the  differentiation 
operator  (matrix). 

For  different  construction  of  the  iVth  degree  polynomial  Q{x),  one  could  have  different 
type  of  boui.  !ary  treatments.  For  examples, 

1.  if  Q{x)  =  0  =  Pc  and  x  =  are  the  (Jauss-Lubatto-(Uiebyshev  points, 

then  we  have  the  ('hebyshev-Legendre  method  ((--L). 

2.  if  Q{x)  =  '  0  =  Or  and  x  =  Xj  are  the  Causs-Lobatto-Chebyshev  points, 

then  we  have  the  ('hebyshev  penalty  method  ((’-P)- 

d.  if  Q{x)  ~  ^ •  <  0  =  Oi  and  x  =  j/j  are  the  (5auss-lA>l)atto- Legendre  points,  then 
we  have  the  Legendre  penalty  method  (L-P). 

Let  d«‘note  e*’‘'  =  e(x^,L,)  and  A/  be  the  time  step  increment,  then  for  j  =  0 . V. 

we  woi'.ld  advance  the  system  of  ODE  (5.2)  in  time  by  the  third  order  Heun  Kunge  Kutta 
scheme  that  has  the  following  form  ; 

For  j  =0.1,...,  .V,  and  =  f/{0) 


•  i 


=  +  y  (O-’!"'  -  rQUM"'  - 

=  + x*""'"  ■  '■‘^•■"'>1''!"  -  "''“I  - 

-  nil,)  -  “?'(/.)  - 


•  { 


wliere  atid  g"(tn)  are  the  derivative-:  of  the  titr.e-dependeiit  boundary  eonditions  in 

time  at  t  = 

It  has  been  observed  before  that  if  one  imposes  boundary  eondition  at  each  intermediate 
stages  of  tlie  Runge  Kutta  scheme,  a  larger  time-step  ((!FL  number)  can  be  used.  Otherwise, 
(!FL  number  has  to  be  reduced  by  as  much  as  four  time  for  stability.  In  this  study,  we  define 
M.  =  CFLjNK 

The  traditional  way  of  the  imposing  exact  boundary  condition  at  or  =  1  can  be  described 
as  following  : 


=  //(a0) 

.J" 

— 

o'’" +  T 

vi" 

= 

At 

+  y) 

of’ 

= 

(n)  .  '^At  (1) 

v]  '-f  -y-Ov; 

^0  ■ 

= 

oiLi  +  *  .j  ) 

rf’+rl'’+' 

= 

£?(<„  +  At) 

(5.4) 


3^ 
4  ^ 


However,  as  shown  in  Table  I  that  this  procerlure  would  lead  to  reduction  of  accuracy  in 
time  as  N  increases. 

Table  1 

Li  Error  and  order  of  accuracy  for  (5.4)  with  k  =  I 


N 

Error 

Rate 

Error 

Rate 

Errttr 

Rate 

16 

0.82E-03 

O.lOE-03 

0.29E-05 

32 

0.15E-04 

2.89 

0.18E-05 

2.91 

0.28E-07 

3.35 

64 

0.42E-06 

2.57 

0.49E-07 

2.61 

0.72E09 

2.64 

128 

0.17E-07 

2.31 

0.19E-08 

2.33 

0.28E-10 

2.34 

(’FL 

8 

4 

1  _ 

Hence,  the  above  pri>c<Hlure  is  modifietl  a.s  following  (Detaihsl  tliscu.ssion  anti  analysis 
will  appear  in  a  ftjture  paper)  ; 

Ft»r  j  =  0. 1, ....  .V,  and  ej***  = 


^  3  ^ 
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I 


w 


»o'’  =  s(<h)+ Y!7'(<n) 

♦i 

+  (5.5) 

-S”"  =  a(<«+.) 

We  shall  denote  this  procedure  as  (XBO).  Table  I  indicated  that  the  order  of  time  accuracy 
for  this  procedure  is  third  order  for  all  N, 

Table  II 

L-i  Error  and  order  of  accuracy  for  (5.5)  with  k  =  1 


N 

Error 

Rate 

Error 

Rate 

Error 

16 

0.77E-03 

0.98E-04 

0.28E-05 

32 

0.12E-04 

3.00 

0.15E-05 

3.00 

0.24E-07 

64 

0.19E-06 

3.00 

0.24E-07 

3.00 

0.37E-09 

128 

0.30E-08 

3.00 

0.37E-09 

3.00 

0.58E-11 

CFL 

I  8 

4 

1 

Next,  using  the  (’-L  method,  one  get  the  La  error  and  the  order  of  accuracy  as  listed  in 
Table  II. 

Table  111 

Li  Error  and  order  of  accuracy  for  (’-L  method  with  k  =  l,r  =  iuo 


N 

Error  Rate 

Error  Rate 

Error  Rate 

16 

0.47E-03 

0.60E-04 

0.28E-05 

32 

0.74E-05  2.99 

0.93E-06  3.01 

0.15E-07  .3.81 

64 

0.r2E-06  3.00 

0.15E-07  3.00 

0.23E-09  3,00 

128 

0.18E-08  3.00 

0.23E-09  3.00 

0.36E-11  2.99 

CFL 

_ ^ 

4 

1 

lablejy. 

Li  Error  of  (’-L  method  for  different  choices  of  r  =  2u;oo  with  k  =  l,('FL  —  ! 


N 

0  =  8 

0  =  2 

o  =  1 

o  =0,9 

o  =  0..5 

16 

0.31  E-05 

0.2KE-05 

0.65E-05 

0.83E-05 

0..32E-0.3 

32 

0.1.5E-07 

0.1.5E-07 

0.15E-07 

0.1.5E-07 

O.IHE-01 

64 

0.23E-09 

0,23E-09 

0.23E-09 

0.23E-06 

unstable 

128 

0..36E-11 

0..36E-11 

0..36E-li 

unstable 

unstable 

•  i 


From  Table  IV,  we  can  see  that  for  r  <  *2u>o,  (^L  becomes  unstable  while  for  r  >  2u;o, 
the  convergent  of  the  scheme  confirms  the  theoretical  prediction. 


Case  2:  Nonlinear  scalar  PDE 

Consider  the  scalar  nonlinear  initial  boundary  value  hyperbolic  equation 

IJf  =  ~2irkcos{2irk{x  +  f)){l  +  8in(27rA:(j:  +  f)))  (5.6) 

<>0 

with  initial  condition 

f/(x,0)  =  2  +  sin(25rl;x) 
and  boundary  condition  at  x  =  1 

=  g{t)  =  2  +  sin(25rA:(  1  +  t)) . 

This  PDE  has  an  exact  solution  given  as  f/(x,f)  =  2  +  sin(2irAr(x  +  f)). 

Table  V 

Li  Error  of  (/'-L  method  for  different  choices  of  t  =  2u;ooi  with  k  =  \,cfl  =  \ 


N 

«  =  8 

rt  =  4 

rt  =  3 

II 

Exact  BC 

16 

0.86E-02 

O.lOE-01 

0.18E-01 

0.27E-01 

0.72E-02 

32 

0.40E-07 

0.40E-07 

0.40E-07 

O.llE-04 

0.39E-07 

64 

0.68E-09 

0.68E-09 

0.68E-09 

unstable 

0.67E-09 

128 

O.llE-10 

O.IIE-IO 

O.llE-lO 

unstable 

O.llE-10 

Different  values  of  k  are  also  tested,  similar  results  are  obtained. 
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